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We present a numerical method for the investigation of non-ergodic effects in the Coulomb glass. 
For that, an almost complete set of low-energy many-particle states is obtained by a new algorithm. 
The dynamics of the sample is mapped to the graph formed by the relevant transitions between 
these states, that means by transitions with rates larger than the inverse of the duration of the 
measurement. The formation of isolated clusters in the graph indicates non-ergodicity. We analyze 
the connectivity of this graph in dependence on temperature, duration of measurement, degree of 
disorder, and dimensionality, studying how non-ergodicity is reflected in the specific heat. 
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I. INTRODUCTION 

Disordered systems of interacting localized particles have been extensively studied for over two decades. A char- 
acteristic feature of these systems is a complex many valley structure of the energy landscape of the state space [EJ . 
Therefore, at sufficiently low temperatures, the system cannot be considered to be in thermodynamic equilibrium: 
^vq Gibbs ensemble theory of statistical mechanics, which is based on the equivalence of time and ensemble averages, is 
\Q . not applicable. Thus non-ergodic effects are important. 

The Coulomb glass |2|J^] is a prominent example of such disordered systems. In heavily doped crystalline semicon- 
ductors, amorphous semiconductor-metal alloys, and granular metals, it plays an important role as a semiclassical 
model for systems of localized states. The dynamical behavior of the Coulomb glass has been studied by several 
groups: Schreiber et al. PH7I, as well as Perez-Garrido et al. [||,[| determined numerically the transition probabilities 
between low-energy many-particle states, and studied the eigenvalues of the transition probability matrix. The former 
j" 1 . group directly diagonalized this matrix, whereas the latter developed a renormalization method to eliminate the tran- 
sitions with large rates, what considerably simplifies the diagonalization. A broad distribution of relaxation times over 
several orders of magnitude was found in both cases. It reflects the glassy behavior of this system. Moreover, Wappler 
et al. ]l(| used the damage-spreading algorithm to study the temporal evolution of the system, and found evidence 
for a dynamical phase transition. Yu studied the time development of the Coulomb gap considering a self-consistent 
equation for the density of states jll| . She too observed that very long time scales are involved. 

The main aim of this work is to study numerically non-ergodic effects in the Coulomb glass analyzing the transitions 
between many-particle states. We apply this procedure to the investigation of such effects in the specific heat. The 
paper is organized as follows: Section II introduces the Coulomb-glass model. Section III describes how the low- 
energy many-particle states are obtained numerically. In Sec. IV, we calculate the transition probabilities between 
these states, and map the dynamical behavior of the Coulomb-glass sample to a graph. The nodes of the graph 
represent the many-particle states, and the edges the relevant transitions between them. Analyzing the structure of 
this graph, we determine the value of a physical observable in dependence on the duration of its measurement. In Sec. 
V, we use this method for the investigation of the non-ergodic effects in the specific heat: We study the influence of 
temperature, duration of measurement, disorder, and dimensionality. Finally, in Sec. VI we extract some conclusions. 



II. MODELS 



The classical impurity band (CIB) model is the most realistic model for simulating an impurity band of localized 
states in a lightly doped semiconductor when quantum interference can be neglected P,p|,[l2|, It is applicable if 
the following two conditions are fulfilled: (i) The mean nearest-neighbor distance is considerably larger than the 
localization radius of the wavefunction of an isolated impurity state, (ii) The temperature is so low that both the 
activation to the conduction band / from the valence band, and the formation of doubly charged donors / acceptors 
can be neglected. 
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We consider a d-dimensional sample of an n-type, partially compensated semiconductor with donor concentration 
Nt>, and acceptor concentration Na = KN^. The degree of compensation K can range from to 1. So donors are 
either occupied, that means neutral, or empty, that means positively charged. Acceptors captured an electron each, 
and are negatively charged. The distribution of electrons between the donors is governed by the Hamiltonian 



H= 1 - - i V — ~ — (l) 



The donor occupation number m equals 1 for occupied donors, and for ionized donors. Moreover, n v = — r„| and 
fij = —Tj\, where the random positions of the donors are denoted by and r,-, and those of the acceptors by r v . 
However, to minimize size effects, we impose periodic boundary conditions and use the minimum image convention p3[ . 
That means, in computing |rj — r,-|, we substitute the projection of |r^ — r 3 -| onto each of the coordinate axis, 

with 8 — 1, by the smallest related value in a periodically repeated representation, Min(|;q — Xj \, L— \x^ ~ x j\) 
with L being the edge size of the sample. For numerical reasons, we construct the samples so that the nearest neighbor 
distance exceeds 0.5, a well justified approximation for amorphous semiconductors, but not for crystalline systems. 
In this work, the unit of distance is defined by the donor density p=l, and electron charge, dielectric constant, and 
Boltzmann constant are taken to be unity. 

Within the present study, most of the calculations are performed for a simplified version of this model. Following 
Refs. mm, see also |2|]3| , we consider a partially filled band of elementary charges (particles) localized on a regular 



lattice formed by the 2Vd donor sites. Here, the acceptors are substituted by background charges —K at each of the 
lattice sites, guaranteeing electro-neutrality on average. The disorder is simulated by a random potential e;. Its values 
are uniformly distributed between —B/2 and B/2. Thus, the influence of the randomness of the donor positions is 
ignored, as well as the correlations between the values of the acceptor potential at neighboring donor sites. Moreover, 
the rectangular distribution of the is a simplification neglecting contributions from particular close pairs of donors 
and acceptors. This model is represented by the Hamiltonian 

H = l^ e * n * + 2^ : > ( 2 ) 

■ . - '17 

I l<J J 

where rii £ {0, 1} denotes again the occupation number of site i. As above, is the distance between sites i and j 
according to periodic boundary conditions. The lattice spacing is taken as unit of distance. 

A mixed form of both the models (1) and (2) is obtained in the following way: The sites are positioned at random, 
and the acceptor potential is substituted by a random on-site potential plus the potential of neutralizing charges —K 
at each site ]l6| , p^7[ . This mixed form is more realistic than Eq. (2): it keeps the donor disorder as the Hamiltonian (1), 
but simplifies the disorder contribution from the acceptors. By means of e,-, the influence of the random surroundings 
of host atoms in an amorphous semiconductor can be simulated. This model is in the following referred to as random- 
position- with-random-potential model. 

The relaxation procedures, which we use, alternatively simulate the sample to be isolated, or to be in contact with 
a particle reservoir Jl8[ . The latter means that, instead of H, the grand canonical potential, h = H — ^^^j, is 
minimized. The value of the chemical potential fj, depends on K, i.e., it is fixed by the electro-neutrality condition. 
Here, we first obtain /i performing a canonical simulation with reduced accuracy. Then we calculate the set of 
low-energy states, as described below, taking into account particle exchange with the reservoir. 



III. DETERMINATION OF A SET OF LOW-ENERGY STATES 



For studying low-temperature properties treating correlations exactly, we need to know a set S of almost all many- 
particles states in a certain energy interval above the ground state energy. If the occupation numbers are known, 
the energy of a state can easily and directly be determined from Eq. (1) because we are treating a classical system. 
The problems, however, consist in the binomially large number of possible configurations, and in the existence of 
many local minima. Thus it is a complicated task to obtain such a set of low-energy many-particle states. It has been 
approached by several methods: Mochena and Pollak JlSj ] developed an approximative renormalization-like procedure. 
Schreiber and Tenelsen |^,|f[ used the Metropolis algorithm to collect low-lying states which seems to be favorable 
in comparison to the previous method p2[. Mobius and Pollak ^3|, and Perez-Garrido et al. [ 
algorithms. In the first stage, they obtained sets of local minima by means of relaxation. For that, a complete search 
considering rearrangements of the site occupations including up to four sites, and an incomplete search concerning 
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more complex rearrangements, built up of several low-energy one-electron hops (shifts), were performed in Refs. p3| ] 
and |jl7), respectively. In the second stage, both groups completed the table of low-energy states by systematically 
investigating the neighborhood (in the configuration space) of each of these states. This neighborhood is defined via 
the accessibility within only one of the considered rearrangements. 

Here, we find the set S of N low-energy many-particle states by means of a three-stage algorithm; for a short 
preliminary description see Ref. p6| . In the first two stages, we create, and improve a "backbone" of S, formed by 
metastable states, the number of which can be considerably smaller than N. Then, in the third stage, we complete 
S by systematically investigating the neighborhood of the states found. Our procedure, which includes sophisticated 
local search |24| , and thermal cycling p5j , is explained in detail in the following. 

In the first stage, creating the backbone of S, we repeatedly start from states chosen at random, and simulate 
quenching the sample (i.e., a rapid relaxation) by means of a local search procedure: In an iterative process, we search 
the neighborhood of the actual state for states of lower energy, and accept always the first such state found. The 
process stops when no lower neighboring state exists. 

Our local search algorithm ||24| ensures stability with respect to rearrangements concerning one up to four sites. 
Making use of the branch-and-bound idea in order to avoid unnecessary attempts to a large extent, it considers the 
following rearrangements: 

(a) transition of one electron between the sample and a reservoir, 

(b) arbitrary one-electron hops within the sample, 

(c) rearrangements by performing simultaneously an (a) and a (b) transition, and 

(d) arbitrary two-electron hops within the sample. 

To ensure high efficiency of the simulations, the searches in (c) and (d) have to be restricted to a certain number of 
neighbors. For one-, two-, and three-dimensional systems, we consider the first 4, 8, and 26 neighbors, respectively. 
(Though the number of possible rearrangements increases rapidly with the number of neighbors considered, the portion 
of long-range hops among the energy decreasing rearrangements is small due to the decreasing interaction strength.) 

The second stage consists in extending, and improving this set of metastable states by thermal cycling p5| . For 
that, a further set Ai of metastable states is considered. Initially, it equals a subset of S, containing the states of 
lowest energy. We cyclically choose one of the states from Ai, and apply to it a Metropolis process with a certain 
temperature. This process, which is referred to as heating, is terminated, however, after a small number of successful 
steps. Thus a great part of the information on the ground state, gained within the previous cycles, is retained. Then 
we quench the sample by means of the above local search algorithm. If the energy of the final state is lower than the 
energy of the initial state, the latter is substituted in Ai. 

For each of the final states of these cycles, we check whether or not it is already contained in S. If not, and if, 
moreover, the total number of states is smaller than its maximum value, we add it to S. However, if the final state 
is not contained within S, but S has been already "filled" , we substitute the final state for the state in S of highest 
energy, provided the energy of the final state is lower than the energy of the latter. 

Several details of thermal cycling should be mentioned: When starting, the temperature of the Metropolis process 
is chosen to be equal to the total energy gain in quenching a random state divided by the number of sites. In the 
course of the thermal cycling procedure, this temperature is reduced gradually. The cyclic procedure continues until 
its efficiency becomes so low that, within a "reasonable" number of cycles, the set Ai has not been improved further. 
In the heating, we consider only two kinds of rearrangements, namely (a) and (b). For d = 1, 2, and 3, the latter 
are restricted to the first 4, 8, and 26 neighbors, respectively. During this process, we fix the occupation of those 
sites which have the same rii value in all states contained in Ai because these values obviously favor a low energy. 
Moreover, in order to diminish the portion of unsuccessful attempts in the Metropolis process, we further reduce the 
set of rearrangements considered: We tabulate the relevant rearrangements after each temperature change, as well as 
after adding a new state to Ai. As relevance criterion, we use the condition that, for at least one of the states in Ai, 
the related energy change must be smaller than the temperature. In heating, we consider only the excitations stored 
in this table. Each heating is terminated after having performed 50 rearrangements. 

In the third stage, we complete the set S of low-energy states by systematically investigating the neighborhood 
(in the configuration space) of the states found |l7|,^3|: For each state in S, we construct all neighboring states, 
considering the same types of rearrangements as above. However, we perform (c) and (d) searches (restricted also 
here to the first 4, 8, and 26 neighboring sites if d = 1, 2, and 3, respectively) only for those states, which are local 
minima with respect to the (a) and (b) rearrangements. Each new state is added to S if the number of states in iS is 
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smaller than its maximum value. Otherwise, we substitute the highest-energy state in S provided the energy of the 
new state is lower. 

This completion procedure starts with the state which has the highest energy, and proceeds cyclically until all states 
contained in S have been investigated. If the number of states in S is large, the decision whether or not the current 
state has been added to S already previously, is particularly CPU time consuming. This decision is significantly 
accelerated by searching hierarchically in an array ordered according to the energies, where pointer arithmetics is 
used. 

To check our program, we performed a series of tests. The degree of completeness of the spectrum of many-particle 
states within a certain energy interval above the ground state energy was judged in a similar way as in Ref . Q . We 
observed that the completeness and the CPU time required do not only depend on the number N of states included in 
iS, and on the number of local searches performed in the first stage, but to a considerable extent also on the number of 
states in M.. In the numerical experiments described in Sec. V, N has values between 25 000 and 75 000 what ensures 
that the width of the related energy interval exceeds the temperature by at least a factor of 25. We found that it is 
a good choice to adjust the number of local searches in the first stage to N/25, and the number of states in Ai to 
AT/500. Naturally, the kinds of rearrangements considered in the three stages are very important too. 

We devoted special attention to testing the thermal cycling part, which improves the backbone of <S. Its efficiency in 
finding states of particularly low energy is illustrated by Fig. 1. This graph presents the mean energy of the lowest state 
found in dependence on the CPU time required, contrasting results for simulated annealing (same rearrangements 
as in heating), sophisticated multi-start local search (repeatedly quenching states chosen at random by means of 
rearrangements (a) - (d)), and thermal cycling. Fig. 1 shows that, in searching for states of very low energy, thermal 
cycling is not only far superior to simulated annealing, but also to the sophisticated multi-start local search algorithm. 
Thus the incorporation of thermal cycling leads to a considerable efficiency increase in the construction of sets of low- 
energy many-particle states. 

As a check of our numerical approach to obtaining almost complete sets of low-energy many-particle states, we 
verified that the corresponding results for the equilibrium specific heat agree with Refs. [ p3| , p6[ . 



IV. INFLUENCE OF THE DURATION OF MEASUREMENT 



Our main aim is to study the influence of the duration r m of the measurement, during which the state of the 
sample travels randomly through its configuration space, on the expectation value of an observable O. For that, 
we decompose the configuration space into separate regions, isolated from each other on the time scale r m , in other 
words, into clusters of states. Two approximations are basic for this approach: (i) The duration of measurement is 
long enough for establishing thermal equilibrium inside each cluster, (ii) Transitions between different clusters are so 
rare that they can be ignored on the time scale of the measurement. 

First, we calculate the transition probability per unit time from the many-particle state / to the state J. According 
to Ref. ||, it is given by 



Wi^j = w exp 



exp[{E! - Ej)/T\ 



if Ei < Ej 
if Ei > Ej 



(3) 



In this equation, the parameter wq is a constant of the order of the phonon frequency, wq 
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denotes 



the localization radius, and Ej the energy of the state /. The sum term is an abbreviation of the minimized sum 
of the related hopping lengths; it concerns only those sites, the occupation of which is changed in the transition. 
That means, we decompose the many-electron transition into independent one-electron hops. Among all possible such 
decompositions, we choose that for which the sum of the lengths of the one-electron hops takes its minimum value. 

Presuming the occupation probabilities of the states to have their thermodynamic equilibrium values, we obtain 
the rate of transitions from state I to state J. 



Ri^j = Wi^j exp 



-Ei 
T 



(4) 



where Z is the partition function. In thermodynamic equilibrium, Ri^j — Rj^i. Thus the transition time, i.e., the 
inverse of the equilibrium transition rate, is given by 



ti.j = r exp 



Eu 
T 



(5) 
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with Eu = max(_E/, Ej), and To = w^ 1 . 

Now, we map the dynamics of the Coulomb glass sample simulated to a graph, where the nodes represent the 
states, and the edges those transitions between them for which tjj < r m . The nodes which are directly or indirectly 
connected with each other by such edges form clusters. These cluster correspond to regions of the configuration space 
being isolated from each other on the time scale r m . Here we presume equilibrium in determining 77 j. However, this 
principle can easily be generalized to the non-equilibrium case, see below. 

Provided, at the beginning of the measuring process, the sample is in one of the states of the cluster a, we measure 
as value of the observable 0, 



<0> a (T) = £0 7 exp(^) Z~\ 
lea ^ ' 



(6) 



where Oi denotes the value of O for the state /, and Z a the partition function of this cluster. The values of {0) a can 
have a very broad distribution. 

Finally, we turn to the expectation value of O, that is the mean value of repeated measurements, referred to in the 
following as ((0))(T, r m ). The probability to find the sample in one of the states belonging to the cluster a be P a . 
Thus, ((O)) is given by the weighted average of the (O a ): 

((0))(T,r m ) = 2<0> Q P Q . (7) 

a 

In consequence, ((O)) depends via the cluster structure and P a on r m , and also on T. 

Both the cluster structure and P a are influenced by the history of the sample. In our numerical experiments, we 
have simulated two situations: 

(A) As standard case, we presume that the sample has reached thermal equilibrium before the measurements; thus 
Pa = Z a /Z . 

(B) Alternatively, we assume the sample to have been quenched from infinite T to the measuring T within a short 
time interval r q . To emulate this process we quench first to T — 0, and heat then immediately to the measuring 
T: In the beginning, we assign the same probability to all states, Pi = 1/N. Then, we perform the following 
iteration treating the states according to decreasing energy. We re-distribute the weight of the considered state 
to all those states of lower energy to which a transition can happen within r q modifying the set of the Pi 
correspondingly. The related transition times are obtained in a way analogous to the derivation of Eq. (5). That 
means, substituting the equilibrium occupation probabilities by the (iteration cycle dependent) non-equilibrium 
Pi, we get for the transition I — > J 

f 00 if Ei < Ej 

TW 1 ^o 1 exp(2E^/a) Pf 1 X Ej > Ej ' { * > 



Thus the quenching process is related to the modification of parts of the states graph (differing from the 
"equilibrium graph") in each of the iteration steps. At the end of the iteration, only the local minima have 
a finite occupation probability. Finally, we assign to each "equilibrium cluster" the sum of the occupation 
probabilities of the included "non-equilibrium local minima" . In this way, we calculate the non-equilibrium 
values of P a , but we neglect the influence of the preparation mode on nj, and thus on the structure of the 
clusters which are isolated from each other on the time scale r m . 



V. NON-ERGODIC SPECIFIC HEAT 



Applying the methods presented in Sees. Ill and IV to the study of non-ergodic effects in the specific heat, we start 
from the definition of this observable for a cluster a of states, which is in thermal equilibrium: 

Since the specific heat is by itself a quantity relating to an ensemble of states rather than to a single state, we do 
not need brackets here to mark the averaging over the states in the cluster a, in deviation from the general notation 
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(0) a (T) in Eq. (6). After averaging over the set of clusters, the r m dependent value of the specific heat (c)(T, r m ) is 
obtained utilizing Eq. (7). 

To make the influence of r m directly visible, we consider the quotient of the values of the specific heat for finite and 
infinite duration of measurement, respectively: 

9 (T,r m )= < / c) ^' Tm . ) . (10) 
(c)(T, oo) 

This quantity is in the following referred to as non-ergodicity quotient. Studying q(T,r m ) rather than (c)(T, r m ) is 
additionally motivated by numerical reasons: The random fluctuations of q from sample to sample are considerably 
smaller than the fluctuations of (c). Moreover, in order to characterize the fluctuations of c from measurement to 
measurement, we consider its root-mean-square deviation, related to the distribution of the clusters, 

a(T,T m )-^P a { j . (11) 

Investigating the physical properties of macroscopic systems, we have calculated ensemble averages, and have 
compared the results for different sample sizes. The number of samples to be taken into account depends not only on 
the accuracy to be achieved, but also on all the various details of the simulated situation. Generally, we have chosen 
the size of the samples so large that corresponding convergence of the results is ensured. The related parameter values 
are given in the figure captions. 

Several details of our simulations have to be mentioned: The localization radius a equals always 0.2 |27|| , Decreasing 
its value corresponds, crudely speaking, to stretching the time scale, and thus enhances the non-ergodic effects. Since 
a great part of the literature concerns the simplified Coulomb glass model according to Eq. (2) (sites arranged on 
a regular lattice), we, too, consider this situation in most of the simulations. Thus, if not stated otherwise, our 
numerical experiments have been performed for this model. The results presented here have been obtained by means 
of a canonical procedure: In calculating (c), we take into account only the low-energy states with total charge equal 
to KNd. When the grand canonical procedure is used for this aim, finite-size effects are more important since, due 
to differing total charge, some clusters remain separated from each other in configuration space even if r m — > oo. 

In the following, we discuss our numerical results in detail. First, we consider the non-ergodicity quotient q. As 
examples, Figs. 2a and 2b illustrate for three-dimensional samples with two different degrees of disorder, how finite- 
size effects influence the dependence of q on the duration of measurement r m . They show that finite-size effects are 
particularly important for large r m , where extended clusters of states have to be considered. However, if the sample 
size Ad exceeds a certain value, roughly 500 for the cases considered here, q is almost independent of Ad. Certainly, 
this limit depends on the concrete situation considered, i.e., on T, B, and d, and also on r m . 

Moreover, Fig. 2 leads to an important conclusion: The value of r m where q reaches (almost) 1 can be interpreted 
as relaxation time of the specific heat. It exceeds typical experimental times by several orders of magnitude. This 
indicates glassy behavior, where non-ergodic effects are fundamental. 

In Fig. 3, we show g(r m ) and its fluctuation a defined by Eq. (11). These data illustrate that the measured values 
of c are distributed within a wide range around the mean value, what originates from the varying properties of 
the separated clusters of states. This is a manifestation of the non-ergodic effects. However, for drawing further 
conclusions, a detailed study of the size dependence of a would be needed. 

Figure 4 displays the temperature dependence of the specific heat for two values of r m , as well as for thermal 
equilibrium (infinite r m ). Moreover, this figure compares the changes caused by variation of T with the influence of 

T m . 

Now we study the non-ergodicity quotient q in more detail. Fig. 5 illustrates that the strength B of the disorder 
has only a weak influence within the considered parameter range. Only for B = 1, a small systematic shift of g(r m ) 
towards higher r m is detectable in comparison to the curves for stronger disorder. If r m is so large that q reaches 
almost 1, decreasing the disorder enhances the relevance of long-time correlations. 

The influence of the temperature T on q is considered in Fig. 6. For r m /ro exceeding 10 10 , we observe the expected 
behavior, namely that decreasing T is related to the increasing influence of non-ergodicity effects. However, we found 
a point, T m /ro ~ 10 10 , where q is almost independent of T within the T region studied. For smaller r m , there is 
even some increase of q with decreasing T. Presumably, the reason of this feature is the following. Decreasing T 
influences the cluster structure in two ways. On the one hand, it causes a division of clusters due to the decrease of 
transition rates, and thus an increase of the cluster number. On the other hand, however, the number of really relevant 
clusters decreases due to the decreasing occupation probability of excited states. We suppose that the second effect is 



G 



dominating for small r m . There is some analogy to uncorrelated hopping conductivity, which increases exponentially 
with T in the limit u> — > 0, but decreases as l/T at high frequency p3). A similar behavior has also been found in 
the frequency dependence of the dielectric susceptibility p3] . 

The question to which extent the obtained results are model dependent is answered in Fig. 7. For the chosen 
disorder strength, the results for the CIB model and the random-position-with-random-potential model are almost 
identical. For the lattice model, however, the g(r m ) curve is a bit steeper, but the relaxation time of the specific heat 
is roughly the same as for the other models. Thus, in the region of low q, its q(r m ) curve is shifted towards higher r m . 
This effect probably originates from the missing of "easy" hops between closely neighboring sites. Thus the influence 
of the spatial disorder of the donor sites is clearly stronger than the influence of the details of the random on-site 
potential e$. (In the CIB model the acceptors create a potential which is almost Gaussian with a width of 2.5.) 

We have also studied the non-ergodicity quotient q for one- and two-dimensional samples. The results, presented 
in Fig. 8, are similar to the above findings for the three-dimensional case. However, there is a clear trend of non- 
ergodicity effects becoming more and more pronounced with decreasing dimension in the region of large r m . Thus the 
reduction of the dimension seems to be related to a stronger "localization" of the particles. 

Finally, we consider the question how the sample preparation influences q(r m ). Fig. 9 shows results for the lattice 
model. Here, we compare samples being in thermal equilibrium with samples prepared by quenching. The latter curve 
is shifted towards larger r m . However, this effect decreases with increasing B, compare with Fig. 2 in Jl^] . Thus it is 
natural that the preparation conditions have a weaker influence for the other two Coulomb glass models considered. 

VI. CONCLUSIONS 

In the previous sections, we have presented a numerical algorithm for studying the non-ergodic effects in the 
Coulomb glass, which are very important at low temperatures. This method considers many-particle states, and takes 
into account correlations completely. Here, it has been used for investigating how the mean value of the specific heat 
(c) depends on the measuring time. The main results of our numerical experiments are summarized in the following 
points: 

(i) The non-ergodicity quotient q(T,T m ) = (c(T, r m ))/(c(T, oo)} strongly depends on the measuring time r m . It 
vanishes as r m — > 0, and, in the cases studied here, it approaches 1 only for values of r m which exceed realistic 
measuring times by orders of magnitude. 

(ii) For large r m , q decreases with decreasing T, whereas an increase of q with decreasing T is observed in the region 
of small r m . 

(iii) Spatial disorder has a larger influence on the non-ergodic effects than the strength of the on-site random potential 
in the lattice model. 

(iv) The importance of non-ergodic effects increases with decreasing dimensionality. 

(v) Preparing the sample by a quench causes an increase of non-ergodic effects in the case of weak disorder. 

Finally, we would like to stress one technical aspect of the simulations. We have improved our previously developed 
algorithm for the construction of almost complete sets of low-energy many-particle states pgfl by adding a thermal 
cycling |25}| part. This method, originally applied to the traveling salesman problem, combines ideas from Monte 
Carlo and local search algorithms. Here, we have used it when searching for deep local minima in the configuration 
space of the Coulomb glass. Also in this case, thermal cycling has been proved to be highly efficient. 
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FIG. 1. Relation between CPU time, tcpu (in seconds), and mean deviation of the energy of the lowest state found from 
the ground state energy, SE = Em Can — -E gr ound state, for one realization of the three-dimensional lattice model with B = 1 and 
JVd = 1000. X: simulated annealing; A: sophisticated multi-start local search; •: thermal cycling. The values of tcpu relate 
to one 180 MHz PA8000 processor of an HP K460. For simulated annealing and multi-start local search, averages were taken 
from 20 runs, for thermal cycling from 100 runs. In thermal cycling, the ground state was always found within 500 seconds. 
Error bars (la region) are presented if they are larger than the symbols. The interpolating lines are guides to the eye only. 

FIG. 2. Influence of the sample size No on the dependence of the non-ergodicity quotient q on the duration of the 
measurement r m for the three-dimensional lattice model, (a) B = 1, T = 0.006; (b) B — 4, T = 0.008. o, +, and • correspond 
to No = 216, 512, and 1000, respectively. The ensemble averaging took into account 200 samples in each case, except for the 
two curves for No = 512 and 1000 in (a), where we considered 100 and 50 samples, respectively. The error bars are smaller 
than the symbol size. 

FIG. 3. Fluctuation of the contributions to q from the individual clusters, represented by error bars displaying a defined by 
Eq. (11). B = 1, T — 0.008, and No = 512. Ensemble averaging took into account 100 samples. 

FIG. 4. Temperature dependence of the mean specific heat for thermal equilibrium (o), and for two finite measuring times: 
Tm/ro = 10 15 (A) and r m /r = 10 12 (x). The error bars represent the width of the distribution of c a . d = 3, B = 1 and 
No = 512. Ensemble averaging took into account 100 samples. 

FIG. 5. Influence of the disorder strength B on g(r m ) for the lattice model: •: B = 1; +: B — 2; o: B — 4. For all curves, 
d = 3, No = 512, T = 0.008. Ensemble averaging took into account 200 samples for B — 2 and B = 4, and 100 samples for 
B = 1. 



FIG. 6. Influence of the temperature T on q(r m ) for the lattice model: •: T = 0.006; A: T = 0.008; +: T = 0.010; o: 
T — 0.012. For all curves: d = 3, B = 1, and A^d = 512; ensemble averaging took into account 100 samples. 

FIG. 7. Comparison of the q(r m ) obtained for the three different versions of the Coulomb glass model introduced in Sec. 
II: +, o, and • denote lattice, CIB, and random-position-with-random-potential models, respectively. For the first and the last 
model, B — 2. In all cases, d = 3, T = 0.008, and No = 512; ensemble averaging took into account 200 samples. 

FIG. 8. Comparison of the q(j m ) obtained for lattice models of different dimensions: o: d = 1, A^d = 50; +: d = 2, 
./Vd = 484; •: d = 3, ./Vd = 512. In all cases, B = 1, and T = 0.01. Ensemble averaging took into account 200 samples for d = 1 
as well as for d = 2, and 100 samples for d = 3. 

FIG. 9. Influence of sample preparation on q{r m ) for the three-dimensional lattice model: o: sample in thermal equilibrium; 
• : sample quenched within r q = 10 13 To ~ 1 s. For both curves: B = 1, T = 0.01, and No = 512; ensemble averaging took into 
account 100 samples. 
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